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ABSTRACT 

We perforin a detailed modelling of the post-outburst surface emission of the low 
magnetic field magnetar SGR 0418-1-5729. The dipolar magnetic field of this source, 
B = 6x10^^ G estimated from its spin-down rate, is in the observed range of magnetic 
fields for normal pulsars. The source is further characterized by a high pulse frac¬ 
tion and a single-peak profile. Using synthetic temperature distribution profiles, and 
fully accounting for the general-relativistic effects of light deflection and gravitational 
redshift, we generate synthetic X-ray spectra and pulse profiles that we fit to the 
observations. We find that asymmetric and symmetric surface temperature distribu¬ 
tions can reproduce equally well the observed pulse profiles and spectra of SGR 0418. 
Nonetheless, the modelling allows us to place constraints on the system geometry (i.e. 
the angles ill and ^ that the rotation axis makes with the line of sight and the dipo¬ 
lar axis, respectively), as well as on the spot size and temperature contrast on the 
neutron star surface. After performing an analysis iterating between the pulse profile 
and spectra, as done in similar previous works, we further employed, for the first time 
in this context, a Markov-Ghain Monte-Garlo approach to extract constraints on the 
model parameters from the pulse profiles and spectra, simultaneously. We find that, 
to reproduce the observed spectrum and flux modulation: (a) the angles must be re¬ 
stricted to 65° ^ ■*/' + C ^ 125° or 235° ^ '0 + C ^ 295°; (b) the temperature contrast 
between the poles and the equator must be at least a factor of ~ 6, and (c) the size 
of the hottest region ranges between 0.2-0.7 km (including uncertainties on the source 
distance). Last, we interpret our findings within the context of internal and external 
heating models. 

Key words: pulsars: general - stars: magnetar - stars: magnetic field X-rays: indi¬ 
vidual: SGR 0418-f5729. 


1 INTRODUCTION 

Isolated neutron stars (NSs) are characterized by a bewil¬ 
dering variety of astrophysical manifestations. Among those, 
particularly intriguing is a class of sources characterized by 
long periods (P ~ 2 — 11 sec) and high quiescent X-ray lu¬ 
minosities {Lx ~ 10^® — 10^® erg s“^), generally larger than 
their entire reservoir of rotational energy | |Mereghetti|200^ . 
These sources, historically classihed as Anomalous X-ray 
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Pulsars (AXPs) and Soft Gamma-Ray Repeaters (SGRs), 
often display stochastic bursts of X-rays, releasing energies 
~ 10®® — 10^® erg in timescales of seconds or less, and spo¬ 
radic though very energetic y-ray flares, with typical en¬ 
ergetics ~ lO'^"^ — lO"'® erg. Furthermore, AXPs and SGRs 
also show long-term outbursts, where their X-ray emission 
increase up to several orders of magnitudes in days/weeks, 
and decays on timescales of years (Rea fc Esposito|2011|. 


The most successful model to explain both the high 
quiescent X-ray luminosities, as well as the X-ray bursts, 
giant y-ray flares, and long-term outbursts, is the magne- 
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tar model ( Thompson &: Duncan|[l995 1996), according to 
which SGRs and AXPs are NSs endowed with large mag¬ 
netic fields, B ~ 10^'* — 10 ^^ G, resulting for example from 
an active dynamo at birth (Kouveliotou et al.[|l998|Q After 


a magnetar is born, the internal magnetic field is subject 
to a continuous evolution through the processes of Ohmic 
dissipation, ambipolar diffusion, and Hall drift. In the crust, 
magnetic stresses are generally balanced by elastic stresses. 
However, as the internal field evolves, local magnetic stresses 
can occasionally become too strong to be balanced by the 
elastic strength of the crust, which hence breaks, and the 
extra stored magnetic/elastic energy becomes available for 
powering the bursts and flares (Thompson fc Duncan|1995 


1996). Alternatively, the ourbursts could be triggered by in¬ 


stabilities in the external magnetic flux tubes (Beloborodov 
fc Levin|2014 Link|20i4 Lyutikov|2015 ). 


While the magnetar model has been very successful in 
explaining some general features of the triggering mecha¬ 
nism of bursts and flares, the discovery in 2010 of an out¬ 


burst from SGR 0418-1-5729 (SGR 0418 hereafter, van der 


Horst et al. 2010 Esposito et al.| 2010), a NS with dipo¬ 
lar magnetic field Hdip = (6 ± 2) x 10^^ G (Rea et aL]|2010 


2013), lower than those of most magnetars, clearly showed 


that the overall picture was not complete. Following this 
discovery, two more sources showing magnetar-like activity, 
but with “non-magnetar”-like B fields have then been dis¬ 
covered ( Rea et al.|2012 Scholz et al.|2012 Rea et al.|2014 


|Zhou et al.||2014[). 

During the last few years, a number of investigations 
have been aimed at understanding the physical reasons for 


2009 


2011 



[Aguilera et al.||2008a|b[ |Kaspi||20l6[ |P erna fc Pons I 


Vigano et al. 2013). In particular, a suite of magne- 


tothermal simulations highlighted the importance of a hid¬ 
den toroidal field in determining the observational manifes¬ 
tations of a NS ( Pons fc Perna|2011[ Vigano et al.|20l3 ). A 
NS with a low dipolar component of the magnetic field (as 
inferred from timing measurements) could still display an 
outbursting behaviour and an enhanced quiescent X-ray lu¬ 
minosity if endowed with a much stronger internal toroidal 
held ( [Turolla et al.|[2011[ |Perna fc Pons|[2011[ |Rea et ah] 


2012 ). 


However, the presence of a strong toroidal held remains 
hidden from timing measurements. Nonetheless, this com¬ 
ponent of the held leaves its strong imprint on the surface 


temperature of the NS (e.g., Shabaltas fc Lai 2012 Perna 


et al.|[2013 Geppert fc Vigano||2014 ), which can be probed 
by means of phase-resolved spectral analysis of the quies¬ 
cent X-ray emission. The goal of this work is to perform 
such an analysis on the post-ouburst emission of the low- 
B held source SGR 0418 with the aim of constraining its 
surface temperature distribution and, in turn, gaining some 
insight into the topology of the magnetic held in the NS 
crust and into whether an additional (external) source of 
heating may be needed (e.g., see Bernardini et al)]|2011 for 
a similar analysis on the quiescent emission of the transient 
magnetar XTE J1810-197). 


The structure of this article is as follows: Section 
presents the surface emission model used and details the 
computation of the spectra. Section describes the data 
reduction and the spectral and timing analyses performed. 
The results of the modelling are presented in Section]^ hrst 
using an iterative analysis, then with a simulateneous analy¬ 
sis of the pulse prohle and spectrum. In Section]^ we discuss 
our hndings within the context of magnethotermal models, 
and hnally. Section [^concludes this article with a short dis¬ 
cussion and a summary of the results. 


2 SPECTRAL MODELS FOR SGR 0418+5729 

2.1 Family of temperature profiles for 
SGR 0418+5729 

The timing properties of SGR 0418, as well as its X-ray lu¬ 
minosity, are consistent with those of an evolved magnetar 
which experienced substantial held decay but still retains a 
strong enough internal toroidal held. In particulajJTuToha 


et al. (2011), using the magnetothermal code by Pons et al. 


(2009), found that the quiescent luminosity of the source and 


its timing properties are compatible with those of an old NS 
born with a super-strong magnetic held which underwent 
signihcant decay over a time « 10® yr. More recently, follow¬ 
ing a rehned timing solution ( Rea et ar]|2013 ), the realistic 
ag^ of SGR 0418 was estimated to be ~ 550 kyr with the 
state-of-the-art magnetothermal evolution model of |Vigaii6| 
|et al.| ( |2013 ). The initial strength of the dipolar component 


was estimated to be in the range Rdip ~ 1 — 3 x 10^'^ G 
(Turolla et al. 2011 Rea et al. 2013 a larger value would 


have spun down the pulsar too much during its estimated 
lifetime). However, in order to display a non-negligible out¬ 
burst rate, they argued (using the results from the simula¬ 
tions of Perna fc Pons|[2011| and [Pons fc Perna 2011) that 
the internal toroidal held at birth had to be much larger. 

As discussed in Section independent constraints on 
the crustal magnetic topology can be obtained through the 
imprints of the magnetic held on the surface temperature of 
the star. In the NS crust, the coupled evolution of the tem¬ 
perature and the magnetic held gives rise to an anisotropic 
temperature prohle, with the degree of anisotropy being con¬ 
trolled by the ratio between thermal conductivity along and 


across the held lines (for a complete description, see Vi- 
|gano et al.||2013] ). As the NS ages, the magnetic held in its 
crust evolves under the combined inhuence of the Lorentz 


force (causing the Hah drift, see e.g. 

Goldreich fc Reiseneg- 

ger|1992| [Hoherbach & Rudiger|2002 

2004| [Gumming et al. 

2004| Pons et al.|2007 Vigano et al.| 

2012| [Gourgouliatos & 

Cumming|2014 ) and the Joule effect (responsible for Ohmic 


dissipation). For a held at birth which is predominantly 
poloidal, the symmetry with respect to the equator is main¬ 
tained throughout the evolution. However, the presence of 
strong internal toroidal components can drastically change 
this topology. If the toroidal held is dipolar, the equatorial 
symmetry is broken during the evolution due to the Hall 
term in the induction equation (e.g., [Vigano et al.|[2013[) 


^ An updated compilation of the measured dipolar fields can be 
found here: http://www.physics.mcgill.ca/~pulsar/niagnetar/ 


main.html 


I Olausen fc Kaspi|2014 l. 


^ This is to be compared with the charasteristic age = P/2P = 
35Myr of SGR 0418. 
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Distance from pole (km), for RNs=10km 
0.05 0.1 0.2 0.5 1. 2. 5. 10. 



Figure 1. Synthetic profiles of the temperature (in units of Tpoig) 
on the surface of the NS as a function of the angle from the 
magnetic pole (6 = 0°). A range of profiles is shown here, with 
different values of e and n (see Equation]^: Red, Blue and Green 
curves correspond to e = 0.25, 0.50, 0.75, respectively. The dotted, 
dashed and solid curves correspond to n = 5000, 10000, 20000, 
respectively. 


which leads to a complex field geometry with asymmet¬ 
ric north and south hemispheres. Hence asymmetric tem¬ 
perature profiles are produced, with a degree of anisotropy 
strongly dependent on the initial toroidal field strength. On 
the other hand, a strong internal quadrupolar toroidal field 
will maintain the symmetry between the two hemispheres, 
while increasing the temperature contrast between the hot¬ 
ter and coler regions on the NS surface. A suite of magne- 


tothermal simulations by Perna et al. (20131 and Geppert 


& Vigano (20141 illustrated the effects of a strong inter¬ 


nal toroidal field as the NS ages. Temperature differences 
by more than a factor of two between the two hemispheres 
could be produced in evolved objects, though the specific 
quantitative details depend on some elements of the micro¬ 
physics such as the magnitude of the conductivity and the 
composition of the envelope. Also the mass of the star influ¬ 
ences the rate of cooling, as well as the extent by which the 
field penetrates into the core ( Vigano et al.|2()T3 l. 

For the purpose of this work, starting with a large set 
of magnetothermal simulations by varying all the relevant 
parameters described above and trying a one-by-one fit to 
each of them is impractical, especially in light of the fact 
that the magnetothermal simulations are computationally 
expensive. Therefore, we adopt the following strategy: 


(a) We use the observed properties of the pulsed pro¬ 
file, together with the qualitative predictions for the field 
strength/configuration expected for an evolved magnetar 
with the characteristics of SGR 0418, to produce a family of 
analytical, parameterized temperature profiles, with which 
spectra and pulsed profile are fitted. The fits largely restrict 
the allowed parameter space for the temperature profiles. 

(b) We then discuss our results within the context of mag¬ 
netothermal simulations, as well as possible external heat¬ 
ing, to identify the most likely physical scenario for the pro¬ 
duction of temperature profiles consistent with those derived 
from the fits. 


Single-pulse profiles with high pulse-fractions (PFs) 


such as that resulting from the observations of SGR 0418 
cannot be readily explained by pairs of antipodal hot spots 
on the surface of the rotating NS. In fact, |Beloborodov| 
(20021 showed that two isotropically emitting, infinitesimal 
hot spots on the surface of a NS (of 1.4 solar masses and 
10 km in radius) cannot produce PFs larger than about 12%. 


Poutanen & Beloborodov (20061 generalized the calculation 


to allow for anisotropy in the local emission, and obtained 
PFs of at most a few tens of percent. Similarly, |DeDeo et al.| 
(20011 studied the modulation level of the emitted flux as a 
function of the hot spot size, for two antipodal spots. They 
found that, for a NS of 1.4 solar masses, radius 12.3 km, 
and a temperature contrast of 9 between the spots and the 
rest of the surface, the PF can reach a maximum of ~ 55% 
if the local emission is strongly pencil-beamed, with an in¬ 
tensity oc cos® 5, where <5 is the angle of the emitted pho¬ 
ton with the normal to the surface. However, a more recent 
study was performed by Shabaltas & Lai (20121, using a 
family of parametrized temperature profiles with symmetry 
with respect to the equator. By employing a new method 
to efficiently compute the radiation intensity from different 
patches of a NS surface with arbitrary magnetic fields and 
effective temperatures, they were able to produce synthetic 
X-ray light curves for a variety of geometric configurations. 
In contrast to previous works, they were able to generate 
single-peaked profiles with PFs as high as ~ 60%. They in¬ 
terpreted their findings as the result of diffuse hot spots of 
finite size (with varying temperatures), combined with the 
beaming due to anisotropic photon opacities in magnetic 
fields. We note that, consistently with previous work, their 
single pulsed profiles from symmetric temperature distribu¬ 
tions displayed a plateau around the flux minimum. Using 


the same methodology of Shabaltas & Lai (20121, Storch 


et al. (20141 modelled the high X-ray PFs of the pulsar PSR 


B0943-I-10 by means of two asymmetric hot spots, a smaller 
one with stronger B field, and a larger one with smaller B 
field. Once again, the beaming effect of the atmosphere was 
found to be a crucial ingredient in producing a high modu¬ 
lation level. 


baltas & Lai 

(20121, as well as the pulsar PSR B0943-I-10 

modelled by 

Storch et al. (20141, had relatively well mea- 


sured distances, and in addition, PSR B0943-I-10 had the 
viewing geometry constrained by radio observations. As a re¬ 
sult, the area of the thermally emitting region could be well 
determined. On the other hand, the distance to SGR 0418 is 
rather uncertain. If located in the Perseus arm (which is in 


its direction), then the distance would be about 2kpc (van 
|der Horst et aL]|2010[ ). However, there is no independent 
observational link to such an association. Because of this 
uncertainty, we choose to leave the distance to the source as 
a free parameter, which, as discussed in the following (see 
Section 4.11, will turn out to be correlated with the tem¬ 
perature profile on the NS surface, and especially with the 
size of the spot. In addition, given the characteristics of the 
source, and in particular its symmetric pulse profile (i.e. no 
evidence for a plateau around the minimum), as well as the 
very high PFs (~ 80%, and even consistent with ~ 100% 
above 1.2 keV), we focus our exploration on a family of 
temperature profiles asymmetric with respect to the equa¬ 
tor, with one hemisphere hotter than the other. This choice 
is less restrictive in that it allows a wider range of sizes for 
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the hotter component (and correspondingly a wider range 
of distances). However, we will also explore whether, sim¬ 


ilarly to the cases studied by Shabaltas & Lai (20121 and 
Storch et al. (20141, the viewing/emission geometry com¬ 
patible with the data allows for the presence of a second 
antipodal hot region. 

We begin by exploring the following family of temper¬ 
ature profiles 


T (9) = e Tpoie + (1 - e) Tpoie cos" (^0 


( 1 ) 


where 9 is the azimuthal angle, e represents the fractional 
temperature difference between the two poles, and n is an 
integer power which measures the gradient of the temper¬ 
ature decline between the two poles. Some representative 
profiles from this family are displayed in Figure for dif¬ 
ferent values of the parameters e and n. 

Similarly, the symmetric profiles (i.e., a pair of antipo¬ 
dal spots) are parametrized as 


T {9) = e Tpoie + (1 - e) Tpoie cos" (0) , (2) 


with n constrained to be even. 

We note that there is a certain degree of degeneracy be¬ 
tween the gradient of the temperature profile, the local emis¬ 
sion pattern of the radiation, the compactness Mns / Rns of 
the NS, and the viewing geometry. For the local emission, we 
assume blackbody radiation with an analytical prescription 
for its beaming. We realize that this is an approximation 
in the case that the star has an atmosphere, which pro¬ 
cesses the surface radiation creating a distorted blackbody 
and modifying the emission pattern. The only correct way 
to approach this problem would be to generate a magne¬ 
tized atmosphere model for each patch of the star (different 
B and T). To the best of our knowledge, this has been done 
only with an approximate method to produce realistic pulsed 
profiles coupled with analytical temperature/magnetic pro¬ 
files ( Shabaltas fc Lai|2012 Storch et al.|2014) ). However, it 
is impractical to implement this method while formally fit¬ 
ting data since it would require the computation of a much 
larger number of atmospheric prohles than done in those pa¬ 
pers (we will be minimizing over a wide grid of parameters 
for both the temperature profile and the viewing geometry). 
On the other hand, using a single model atmosphere (i.e., 
a single strength for B and one direction, perpendicular to 
the normal) on the entire surface - which is reasonably easy 
to do in a £t - would be generally incorrect, and especially 
so if the object has a significant non-radial component of 
the external magnetic field. Therefore, given that the de¬ 
tailed magnetic topology of SGR 0418 is not apriori known, 
and the fact that, even with a perfect spectral computation 
there would still remain a degree of degeneracy in the pulsed 
profiles with the compactness ratio of the star (which we as¬ 
sume fixed at some typical value), and the viewing geome¬ 
try, which we constrain by fitting the pulse profile, we adopt 
the simplest approach of assuming local blackbody emission 
with a parametrized form for the beaming which captures 
the ’limb-darkening’ effect of magnetized, light-element at¬ 
mospheres (see also Bogdanov|20 14 for a similar approach). 
However, in order to quantify the dependence of the results 
we obtain on the presence of beaming, we will also repeat 
part of the analysis for isotropic emission. 


2.2 General-relativistic, phase-dependent spectra 


The general relativistic, phase-dependent spectra are calcu- 


lated using the formalism developed by l 

Page|1995 see also 

Pechenick et al. 1983 Pavlov & Zavlin 

2000|. The intense 


gravitational field of the star bends the photon trajectories: 
a photon emitted at an angle <5 with the normal to the NS 
surface will reach a distant observer if generated at an angle 
with respect to the viewing axis, where 



(1 — 2u)u^ 


du. (3) 


In the above equation, x = sin 5, Rs = ^GMjc? is the 
Schwarzschild radius of the star, and R and M are, re¬ 
spectively, its radius and mass. Here we adopt a canoni¬ 
cal mass M — 1.4 M© and radius RNS=10km. These values 
yield a gravitational redshift (1 — Rs/R)^^^ = 0.766. Cor¬ 
respondingly, a distant observer would measure a radius of 
Rrx = R X {I — Rs/R)~^^^ = 13.1km, and a surface tem¬ 
perature of Too = T X (1 — Rs/RY^^■ In the following, we 
will quote £t results for the pole temperature in terms of the 
redshifted values, Tpoie,oo. 

Let Q,{t) be the modulus of the NS angular velocity, and 
let us define the phase angle '){t) = / Q.{t)dt as the azimuthal 
angle subtended by a reference vector n around the axis 
of rotation. As the reference vector, we choose the polar 
axis, in the coordinate system in which T{9 = 0) = Tpoie. 
For a dipolar component of the magnetic field, this would 
correspond to the magnetic axis; for the toroidal component 
of the field, it represents the symmetry axis. 

As the star rotates, the angle between the vector h and 
the line of sight is given by 


a{t) = arccos [cos -i/i cos ^ -I- sin ip sin ^ cos y{t)] , (4) 


where we indicated by ip and ^ the angles that the rotation 
axis respectively makes with the line of sight and n (see Fig¬ 
ure 1 in |Perna fc Gotthelf|2008] for a graphical representation 
of the viewing geometry). Note that there is a degeneracy 
between the two angles %p and i.e., they can be exchanged 
without altering the pulse-profile or the spectra. 

The phase-dependent spectrum measured by the ob¬ 
server as the star rotates is obtained by integrating the local 
emission over the entire observable surface. Accounting for 
the effect of gravitational redshift of the emitted radiation, 
this integral takes the form 

p^2 t) 2 pi p27T 

y y /{T[e„(x),fli„],£;}d<).„da;, (5) 

where D is the distance, Eaa = E{1 — Rs/RY^^ is the en¬ 
ergy measured by a distant observer, the spectral function 
I[T{9v, (pY, E] describes the dimensionless distribution of 
the locally emitted photons (the blackbody function here), 
and {9v,(pv) are the coordinates on the surface relative to 
the line of sight. 

As discussed above, to emulate the effect of a magne¬ 
tized, light-element atmosphere, we assume the local pho¬ 
ton intensity to emerge in a pencil-beaming pattern with 
respect to the normal to the surface, which we model as 
f{S) oc cos^ 5, with a beaming intensity p = 1 as a closer 
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Table 1. XMM-Newton Observations of SGR 0418+5729 


ObsID 

Start 

Usable time 


Time (TT) 

(ksec) 

0693100101 

2012-08-25 14:18:08 

58.7 

0723810101 

2013-08-15 18:06:25 

32.5 

0723810201 

2013-08-17 20:57:17 

36.0 


match to that of realistic, magnetized atmosphere models 
I van Adelsberg fc Lai|2006[ p 

From Equation if: other useful quantities for compar¬ 
ison to observations can then be readily computed. In par¬ 
ticular, the phase-averaged spectrum is given by 

1 

Eave(Eoo) = ^ J F[Eoa,a+)] d'y , (6) 

while the pulse profile in a given (observed) energy band, 

{-£‘1,00 5 ^2,00 } 7 is 

£’(7) = F[E^,ai^)] dEoo. (7) 

'^£^ 1,00 

The PF is defined as 

p-p _ -£max(7) -£min ( 7 ) /ON 

En,ax(7)+-F’min(7) ’ ^ ’ 

where the phases corresponding to the maximum and min¬ 
imum of the flux will generally vary depending on the tem¬ 
perature distribution on the NS surface. 


3 DATA REDUCTION AND DESCRIPTION 
OF ANALYSES 

3.1 XMM-Newton observations 


We use in this work two recent observations of SGR 0418 
with XMM-Newton, acquired in August 2012 and August 
2013 (see Table [^, when the source appears to have ap¬ 
proached quiescence with a stable flux. All observations were 
acquired with the EPIC-pn camera ( Striider et al.|2001 1 in 
Large Window mode, i.e., with a time frame of 48 ms, and 
with the EPIC-MOS cameras ( [Turner et al.|[^Oip in Small 
Window mode, with a time resolution of 300 ms. For the 
timing analysis, data from both pn and MOS cameras are 
used, while only the pn spectrum is used for the spectral 
analysi^ 

Standard reduction procedures are applied, using 
epchain in the XMMSAS vl3.5 package, together with the 
latest calibration files. Photon events times of arrival are 


^ These authors particularly discussed the important effect of 
vacuum polarization on the emergent radiation pattern. They 
showed how, without the inclusion of this effect, the emitted radi¬ 
ation exhibits a characteristic beaming pattern, with a thin pencil 
shape at low emission angles and a broad fan at large emission 
angles. Inclusion of vacuum polarization tends to reduce the gap 
and lead to a featureless, broad pencil beaming pattern (especially 
so for stronger fields). 

^ The uncertainties due to the cross-calibrations of the pn and 
MOS effective areas compensate for the gain in signal-to-noise 
ratio when adding the MOS spectra (see Read et al.|2014| |. 


8 

£ 
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Figure 2. Phase averaged spectral fit with the model described 
in Section]^ (see Equation]^, with n = 12783, e = 0.08 and 
(-0, = (13°, 84°), the most likely values from SectionThe 

spectral shape is mostly driven by the pole temperature Tpoie ^i^d 
the absorption column density while the flux depends on n 
and e with weaker dependence on -i/; and The fit is statistically 
acceptable with y^/dof (prob.) = 1.28/34 (0.13). 


corrected to the Solar System barycenter, using barycen, 
before the phases of all events are calculated using the 
best-fit X-ray timing solution reported in Rea et al. (20131: 
P = 9.07838822 sec, P = 4x10"^® sec s"^ on MJD 54993.0. 
The data were checked for background flares, which were 
removed to limit contamination. The resulting usable expo¬ 
sure times for each observation are listed in Table Finally, 
events are filtered in the 0.3-10.0 keV range with the PAT¬ 
TERN <4 and FLAG = 0 restrictions. The observed flux of 
both observations is consistent with being constant, confirm¬ 
ing that SGR 0418’s luminosity and surface temperatures 
have not significantly changed between observations. 


3.2 Spectral analysis 


Phase-averaged spectra are extracted from the two obser¬ 
vations in the 0.3-10 keV range using 20” circular regions. 
Background spectra are extracted from 90” circular regions 
devoid of X-ray sources. The response matrices were gen¬ 
erated using the tools rmfgen and arfgen. The extracted 
pn spectra are then combined into a single spectra using 
epicspeccombine, after verification that the flux was con¬ 
sistent with being constant. Finally, events are binned for 
the purpose of the phase-averaged spectral analysis with a 
minimum of 20 counts per bins. For the spectral analysis per¬ 
formed in XSPEC V. 12.8 ( |Arnaud|1996p , 3% systematics are 
added in each spectral bin to account for uncertainties in the 
absolute flux calibration of the instrument ( Guainizzi|2014 1. 
The spectral model used, described in Section 2| is modu- 


lated by X-ray absorption using the wabs model I Morrison fc 
McGammonjjigs'S I. The overall normalization factor of the 


fit, convoluted with the uncertain distance to the source, is 
left free to vary. All errors from the spectral analysis are at 
the 90% confidence level. The spectrum and best-fit model 
are shown in Figure and the results are presented in Sec¬ 
tion |4| 
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3.3 Pulse Profile Analysis 

The phases resulting from the folding at the timing solution 
given above are used to generate the pulse profile. Further¬ 
more, because of the low count rates in the observations, 
we limit the analysis to 10 phase bins. The errors in the 
number of counts in each bin are Poisson errors. The av¬ 
erage background number of counts is subtracted in each 
phase bin. For the background subtraction, the background 
regions used have the same size as the SGR 0418 source 
regions, and are devoid of any other detected X-ray source. 

The PF is calculated according to Equation (|^, where 
Fmax( 7 ) and ^ 111111 ( 7 ) are the majcimum and minimum X-ray 
fluxes measured in the phase bins. The PF for the cumula¬ 
tive 0.3-10 keV band is PF = 0.78 ± 0.09. While first us¬ 
ing the single 0.3-10 keV band for the computation of the 
pulsed profile and modelling, we then also considered split¬ 
ting the full energy range into smaller bands, since this car¬ 
ries a higher constraining power for the underlying model. 
However, we found that splitting in more than two bands 
would result in a too low count rate per band. Hence we 
limited the split to only two bands, 0.3-1.2keV and 1.2- 
10.0 keV, which were chosen to have approximately simi¬ 
lar count rates. We found the PF in the lower energy band 
to be PFb.3-i.2keV = 0.62 ± 0.10, while the higher band 
PP 1 . 2 - 10 . 0 keV was consistent with 1.0. As expected, the PF 
in the high-energy band is higher than in the low-energy 
band. 


4 CONSTRAINING THE SURFACE 

TEMPERATURE THROUGH SPECTRAL 
AND PULSE PROFILE MODELLING 

4.1 Iterative fitting of spectra and pulse profiles 

Constraining the temperature profile requires a coupled 
spectral and timing analysis. The pulsed profiles are most 
sensitive to the run of temperature with angle on the NS sur¬ 
face, while the phase-averaged spectra are most sensitive to 
the overall flux normalization (reflected in Tpoie and on the 
size of the spot, for a fixed NS radius), and to the amount 
of absorption, quantified by the column density of hydro¬ 
gen Ah, noted Ah ,22 hereafter when expressed in units of 
10^^ atoms cm“^. Once Tpoie and Ah are measured from the 
spectra, fitting the synthetic pulse profiles (computed via 
Equation to the observed pulse profiles allows to con¬ 
strain the system geometry and the temperature profile. 
Here, we first perform the spectral and timing analysis iter¬ 
atively rather than simultaneously. We will show that this 
is a rather good approximation indeed (see also Bernardini 


et al. 2011). However, Section 4.2 will present and validate 
a simultaneous analysis using a Markov-Chain Monte-Carlo 
approach. 

Using the spectral model described in Section to¬ 
gether with the family of temperature profiles defined by 
Eq. §, we first obtain an estimate of Tpoie and Ah from 
a small selection of surface temperature distributions, for 
a viewing geometry = (90°, 90°). In the first step of 

the iteration the viewing geometry is not constrained yet, 
and the choice of an orthogonal rotator is made as a start¬ 
ing point given the high PF of the source. We And that for 
a fixed distance (e.g. 2 kpc), the spectra strictly constrain 



£ 


Figure 3. Pulsed fractions (PFs) obtained as a function of the pa¬ 
rameters n and e, with the system geometry (i/i, 5) = (90°, 90°), 
which maximizes the PF. Contour lines indicate PFs ranging from 
0.1 to 0.9, with 0.1 increments; the red line indicates the 0.8 PF 
contour. Note that while the low range of n values can produce 
PF such as those observed for SGR 0418, it fails to accomodate 
reasonable distances of this source (see discussion in Section [4^. 


the values of e and n such that the model produces the flux 
observed for this source. This is because these parameters 
relate to the fraction of the NS surface which dominates the 
emission (see Figure]^, and hence to the fit normalization. 
In other words, the shape of the observed spectra constrains 
the average emission temperature and absorption, while the 
flux normalization restricts the size of the spot and the tem¬ 
perature constrast at the surface. We note, however, that the 
distance is estimated solely based on the association of the 
magnetar with the Perseus arm of the Galaxy. This assump¬ 
tion stems from the position of SGR 0418 in the direction of 
the Perseus arm. Should this association be fortuitous, the 
distance to SGR 0418 could be mis-estimated to lower or 
higher values. 

Such uncertainties in the distance translate into less 
constrained parameters e and n. As an example, a spectral 
fit with e = 0.1, n = 10,000 (i.e., a spot size of 0.275 km, 
full-width at half-maximum, FHWM, on a 10 km NS) yields 
a distance of 2.13 kpc, while a fit with e = 0.15, n = 5,000 
(spot size of 0.425 km FWHM) yields 4.9 kpc. However, the 
parameters Tpoie and Ah appear relatively stable against 
change in e and n, or in the viewing geometry. We find 
Tpoie ,00 ~ 0.32 keV and Ah, 22 ~ 0.25 for a wide range 
of £ and n that lead to reasonable values of the distance 
(d~ 0.5-4 kpc), i.e., £<0.10-0.15 and n ~ 5000- 20000. 
Therefore we already have a good starting point in this it¬ 
erative process to constrain the parameters of the system. 
Higher values of £ lead to statistically unacceptable fits. 

We note the existence of a secondary minimum, lead¬ 
ing to acceptable fit statistics, in which Ah ,22 ~ 1.2 corre¬ 
sponds to £ ~ 0.4, for the same pole temperature. This could 
be explained by the fact that the larger averaged flux cre¬ 
ated by the less pronounced contrast (larger e) needs to be 
heavily absorbed to fit the data. However, we reject such a 
large value of the absorption, on account of the fitted values 
of Ah from the high signal-to-noise observations obtained 
during the outbursts | Rea et M^|2013 1. 

With the estimated values of pole temperature and ab- 
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Figure 4. Projected 2D map of maximum p-values when the two 
angle ^ and ^ vary, in the case of beamed local emission (pencil 
beam oc cos(5). The solid line represents the probability p = 0.01. 
The best-fit over the whole 4D-space corresponds to a p-value 
p = 0.05. 



Figure 5. Same as in Figure]^ but for the two bands: 0.3-1.2keV 
(red) and 1.2—10.0keV (blue). Only the p = 0.01 contours are 
plotted to allow comparison between the two bands. 


sorption, we can then explore the e-n parameter space for 
pairs capable of accomodating a PF of ~ 0.80. Figure 
shows the PF as a function of e and n for {ip, = (90°, 90°). 
With such angles, a PF ~ 0.80 can be reproduced for a wide 
range of n, while e remains essentially constant at ~ 0.15. 

We note that the high measured PF of the source is 
partly due to the effect of X-ray absorption. More specif¬ 
ically, absorption can artificially increase the intrisic (un¬ 
absorbed) PF since the low-temperature portion of the NS 
surface (emitting lower energy, less pulsed X-ray photons) is 
more affected by absorption than the hotter parts. The ar- 
tifical increase in PF due to absorption depends both on the 
geometry {ip,i) and the temperature gradient on the staij^ 
The emission model of this analysis can be used to estimate 
the intrinsic PF of SGR 0418, i.e., unaffected by absorption. 
Specifically, for n = 10000 and e = 0.15, the intrinsic PF is 
0.36, while it is 0.88 for Nh ,22 = 0.25, hence quite sensitive 
to A^h- 

Next, choosing ip = 90° and ^ = 90° with the tempera¬ 
ture distribution profile defined by e = 0.15 and n = 10000, 
we perform a detailed spectral analysi^in order to further 
refine the value of the pole temperature Tpoie and the X-ray 
absorption A^h- An acceptable fit to the data is obtained 
with Tpoie,oo = O.SSlo QgkeV and Nh ,22 = O.SGtg p® (see 
Figure]^. However, a better fit is obtained with e = 0.10 
and n = 10000 leading to Tpoie,oo = 0.381 q g® keV and 
Nh ,22 = 0.22l°;°? for Xv/dof (prob.) = 1.28/34 (0.13). The 
flux of SGR 0418 is Tx = 2.0 /)q'4 x 10“^^ erg cm“^ s“^ in the 
0.3-3.0keV range. Although these values of n and e slightly 
overpredict the PF (see Figure]^, the best-fit Tpoie and A^h 

® See jPerna et al.| ( |2000[ for an extensive discussion of the effect 
of absorption on the PFs. 

® For comparison, we provide the effective temperature obtained 
when fitting a bbodyrad model: fcTeU oo = 0.32/jg Q 4 keV and 
Fh .22 = O.lS/jg Qg for x?/dof (prob.) = 1.31/34 (0.10). Adding 
a powerlaw spectral component slightly improves the fit (F-test 
probability ~ 3%). 


are consistent with those estimated for different values of e 
and n in the first step of this iterative process. 

With the values of Tpoie,oo = 0.38 keV and Nh ,22 = 0.22 
obtained from the spectral analysis, we next proceed to con¬ 
strain the system geometry and the temperature distribu¬ 
tion that best reproduces the pulsed emission from the sur¬ 
face of SGR 0418. For the analysis that we describe in the 
following, we first consider the pulsed profile in the full en¬ 
ergy band 0.3-10 keV. Using the model described in Sec¬ 
tion a grid of pulse profiles is generated for multiple sys¬ 
tem geometries (varying the angles ip and ^ between 2° and 
180°, in steps of 2°), and multiple temperature distribution 
profiles (by varying n between 5000 and 20000 in steps of 
500, and e between 0.05 and 0.15 in steps of 0.01). Each of 
the synthetic pulse profiles generated for each set of 4 param¬ 
eters {ip, n, e) is then fitted to the observed XMM-Newton 
pulse profile. From the yf values obtained from these fits, 
projected 2-dimensional maps of maximum p-values are ob- 
tainetfl 

Figure]^ shows the maximum p-value maps for the pair 
of angles ip-P,. The fit to the pulsed profile allowed ns to 
mainly constrain the viewing/inclination geometry. The pa¬ 
rameters n and e were not constrained by the pulse profile 
fitting any further than what was imposed by the spectral 
analysis (via the distance requirement). In other words, n 
and e were essentially unconstrained by the pulse profile fit¬ 
ting in the range 5000-20000 and 0.05-0.15, respectively. 

The best-fit resulting from the grid search in the param¬ 
eter space is obtained for the following parameters: t = 0.10, 
n = 7500 for (V’,0 = = (124°,4°), with x^/dof 

The p-value ( [Vogel et al.|201^ [Tendulkar et al.|2015[ see e.g.,), 
also called the null hypothesis probability in XSPEC, is the prob¬ 
ability of finding by chance a x^ as large or larger than the best-fit 
X^ under the assumption that the null hypothesis is true (i.e., that 
the model does not describe the data). The p-value is calculated 
by integrating the x^ probability density function with the num¬ 
ber of degrees of freedom in the data set between the best-fit x^ 
and oo. 












S. Guillot et al. 



Figure 6. Background subtracted pulse profile of SGR 0418 in 
the energy range 0.3-1.2keV (red) and 1.2-10.0keV (blue). The 
dashed line corresponds to the best synthetic pulse profiles ob¬ 
tained from the grid search described in Section |4.1| for each of 
the two energy bands. 



Figure 7. Same as in Figure but in the case of an isotropic 
local emission. 


(prob.) = 1.9/8 (0.05). However, as can be observed in Fig¬ 
ure]^ a wide range of angles is permitted by the pulse pro¬ 
file fitting (p > 0.01), with very different sets of parameters 
n and e. Figure also shows that the angles and ^ are 
constrained to two narrow bands. Note the observed sym¬ 
metry emerges from the symmetry against exchanges be¬ 
tween 'Ip and ^ (see Equation [^. The most conservative con¬ 
straints on the angles can be summarized by the intervals: 
73° ^ 'i/ + C ^ 140° and 220° ^ V' + ^ ^ 287°. 

A consistency check was then performed by using the 
values of the best-fit angles (124°, 4°) and temperature dis¬ 
tribution parameters (n = 7500 and e = 0.10) in a phase- 
averaged spectral fit. These parameters lead to a best-fit 
temperature Tpoie.oo = O.SdlQ pg keV and a hydrogen col¬ 
umn density Nh .22 = O.SOIq q®, for a corresponding dis¬ 
tance of 0.5 kpc. The pole temperature and the hydrogen 
column density remain consistent with the best-fit Tpoie ob¬ 
tained when fixing = (90°,90°), as done at the be¬ 

ginning of this iterative analysis. In addition, for the best- 
fit temperature distribution found above (n = 7500 and 
e = 0.10), the viewing geometry appears to have an effect on 
the distance (i.e. flux normalization). For {ip,^) = (124°, 4°), 
d = 0.5/lQ;^kpc and for = (90°, 20°), d = 2.7to.9 kpc 

(these uncertainties are essentially flux uncertainties). While 
this was not suspected initially, this dependence of the dis¬ 
tance on the viewing geometry is explained by the fact that 
a highly contrasted spot, more or less-often visible during a 
complete phase depending on the viewing angles, will cor¬ 
respondingly generate a larger or smaller flux in the phase 
averaged spectra. 

The grid search described above was then repeated with 
the pulse profiles computed in the two energy bands 0.3- 
1.2 keV and 1.2-10.0 keV. The resulting best fits pulse pro¬ 
files are shown in Figure and the 2D projections of maxi¬ 
mum p-values are displayed in Figure]^ The analysis of the 
pulse profiles in the two bands shows that the higher en¬ 
ergy bands constrain the angles a little more than the low 
energy band and this is simply due to the larger pulse frac¬ 
tion (consistent with 1.0) observed in the 1.2-10.0 keV band 
compared to the 0.3-1.2keV band. But overall, the param¬ 
eter space is mostly insensitive to the energy band chosen, 
given the low signal-to-noise available in these obervations. 


Last, with the purpose of quantifying the effect that 
the uncertain atmospheric beaming has on our results, we 
repeat the analysis assuming local isotropic emission. This 
results in a narrower allowed parameter space for the angles 
Ip or For example, we And that the angles ^ or ^ must lie 
in the range 104° <ip + ^< 145° and 215° < V' + C < 256°. 
Moreover, the contours are less curved in ip-^ space (Fig¬ 
ure]^. These results can be understood since the effect of 
beaming is that of creating a larger modulation for the same 
parameters. Hence, to reproduce the same level of observed 
modulation with a local isotropic emission beaming, smaller 
viewing angles are not allowed, restricting the angle geome¬ 
try to larger values. 

The best-fit pulse profiles in the locally isotropic [xS /dof 
(prob.) = 1.9/8 (0.06)] and pencil-beamed surface emis¬ 
sion are compared in Figure As it can be seen, the pro¬ 
files are almost indistinguishable. This result stems from 
what just discussed above: when the local emission is 
beamed, the fit parameters adjust so that smaller viewing 
angles/temperature contrast are needed to achieve the same 
observed level of modulation. Therefore, a comparison be¬ 
tween Figure 1^ and yields a quantitative estimate of the 
amount by which the uncertain angular distribution of the 
local emission affects the inferred system parameters. 

4.2 Simulateneous fitting of spectra and pulse 
profiles using a Markov-Chain Monte-Carlo 
approach 

As explained in Section |4.1[ the source spectrum is more 
sensitive to some parameters of the model (e.g.: Tpoie or the 
spot size charaterized by n in Equation]^ while the pulse 
profile is more sensitive to variations in the temperature 
contrast at the surface or in the system geometry {ip or ^). 
Nonetheless, all parameters play a role in modelling the sur¬ 
face emission and its phase- and energy- dependence, and 
Section |4.1| showed that it was somewhat difficult to pre¬ 
dict what variations of a given parameter can have on the 
modelled pulse profile or spectrum. Iterating between pulse 
profile and spectral fits proved to be impractical due to the 
number of parameters. Furthermore, a simultaneous analy- 
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Figure 8. Background subtracted pulse profile of SGR 0418, to¬ 
gether with the best fit synthetic pulse profiles obtained from 
the grid search described in Sections |4.1| and |4.3| (dashed blue: 
single spot with pencil-beamed surface emission, dashed red: sin¬ 
gle spot with isotropic surface emission, dashed green: two spots 
with pencil-beamed surface emission). They are computed in the 
0.3-10 keV band. 


sis allows us to have a global understanding on the roles of 
parameters in the model, and to better understand the ex¬ 
isting degeneracies between parameters. For these reasons, 
a Markov-Chain Monte Carlo (MCMC) approach was then 
employed to sample the parameter space while simultane¬ 
ously fitting the spectra and pulse profiles. In addition, the 
method allows to include known priors on some of the pa¬ 
rameters. 

The implementation of the MCMC algorithm used in 
this work, the “Stretch-Move” algorithm (also called “AfRne 
Invariant Ensemble Sample”, Goodman fc Weare|2010 1, was 
previously employed, described and tested extensively in 


other works (e.g., Wang et al. 2013 Guillot et al. 2013 


Guillot fc Rutledge 20141. It made use of the PyXSPEC 
package, the python version of XSPEC ( |Arnaud|[l996| ), for 
the comparison of the synthetic spectra to the phase-average 
spectra. In this algorithm, multiple chains (called “walkers”) 
are sampling the parameter space and each set of proposed 
parameters is accepted or rejected based on the likelihood of 
the model with the proposed sets of parameters. The next 
proposed step of each chain is then chosen from an affine 
invariant distribution along a line between the current point 
of the chain and the current point of a randomly chosen 
chaiijf] 

In this approach, we used a set of seven parameters 
sampled by the MCMC: 

• the two angles of the system geometry 'tp and ^ (con¬ 
strained to -|- ^ < 180°, as justified by the symmetry, in 
Eq. this is to reduce the size of the parameter space ex¬ 
plored, and the full range can be recovered by symmetry). 

• the surface temperature contrast e (Eq. , 

• the parameter n characterizing the spot size (Eq. [^, 

• the column density of hydrogen Ah , 

• the pole temperature Tpaie, 

• the distance d, for which we choose a fiat prior between 


0.5 and 4kpc (see discussion in Section 4.11. 


® The description of the algorithm and its performance are de¬ 
scribed in I 


Goodman &; Weare 


( 20101 . 


The neutron star radius, which controls the star’s compact¬ 
ness and therefore the amount of light-bending, could also 
be sampled in this approach to potentially obtain some con¬ 
straints on its value. However, it was held fixed here due to 
the relatively low signal-to-noise ratio of the data. Higher 
signal-to-noise data could potentially place constraints on 
the radius, although there exists degeneracies with the an¬ 


gle geometry that would not be fully broken (see e.g., Perna 
fc Gotthelf||2008 Bernardini et al.]|2011 (. 


In this MCMC simulation, 50 walkers were run simul¬ 
taneously for a total of 15000 steps each. After removing 
5000 steps of burn-in, an average acceptance rate of 15% 
was obtained. We also checked for convergence by visually 
inspecting the trace of each parameter. The results of this 
simultaneous analysis of the pulse profiles and spectra of 
SGR 0418 with our model are shown in Eigurej^ as the 1- 
dimensional and 2-dimensional posterior distributions of all 
parameters. 

This simultaneous analysis confirms what has been ob¬ 
served in the iterative method of Section [4.1| In particular, 
the constraint on the geometry of the system (see the tp — ^ 
plot of Figure]^ is well reproduced, albeit with a slightly dif¬ 
ferent shape that results from the coupled constraints from 
the spectra (as noted in Section 4.1 the distance and the an¬ 
gle geometry are related). The constraints on the two angles 
are: 65° < + ^ < 125° or 235° < V’ + C < 295°. 

Furthermore, because a conservative prior on the dis¬ 
tance was included, the size of the spot (parametrized via 
n) is not particularly well constrained. We find that the spot 
size is in the range 1.2°-3.8°, equivalent to 0.2 — 0.7 km for 
the 10-km NS considered here. The spot size distribution 
peaks at 0.35 km. We note that allowing for larger distances 
would therefore include larger spot sizes (i.e, smaller values 
of n) in the posterior distributions. Since this is driven by 
the spectra normalization, it would not affect significantly 
the posterior distributions of Tpoie or Ah. However, a conse¬ 
quence would be a broadening of the “banana”-shape con¬ 
tours of the 2D-distribution of tp-L since a larger spot would 
reproduce the pulse profiles for larger values of the angles pj 
and Note also that because of the symmetry arising from 
the system geometry, the angles in the MCMG run were con¬ 
strained to Ip ^ < 180°, and therefore, the ip — ^ plot only 
shows one of the “banana” shape contours. 

The temperature contrast e is constrained to similar val¬ 
ues as those obtained in the iterative analysis, i.e. e < 0.15. 
However, it is important to note that while allowed here, 
e —>■ 0.0 represents an unphysical description of the temper¬ 
ature distribution at the surface of the NS since the cold 
part cannot have a zero temperature. Furthermore, given 
the posterior 2D distributions of Figure one can readily 
see that excluding the lower range of e, say 0.00-0.05, would 
not significantly affect the other parameters due to the mild 
correlation between e in this range and the other parameters. 

The posterior distribution of the pole temperature is 
Tpoie,(X = 0.36 ± 0.05 keV (90% c.l.) while the hydrogen col¬ 
umn density is Nh ,22 = 0.251q;J 8 (90% c.l.). Note that we 
initially observed the bimodal distribution of Ah and there¬ 
fore e (leading to the secondary yf minimum in the spectral 
fit described in Section 4.11. We then applied a prior on Ah 
to exclude values Ah ,22 > 0.8 (see Section 4.11. The dis¬ 
tance is slightly more constrained than the initially prior 
range allowed. We find a distribution between 1.5 kpc and 
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Figure 9. 2-dimensional posterior distributions for each pair of parameters and 1-dimensional posterior distributions for each parameter 
resulting from the simultaneous fits of the pulse profiles and spectra via the MCMC approach described in Section |4.2| The color scales 
of the 2-dimensional posterior distributions represent the density of points (i.e., blue means zero-density), compared to other figures of 
this article where the color scales represent the p-values. The solid and dashed contours are enclosing 68% and 95% of the accepted 
points. This figure was created with the Mathematica package LevelSchemes ||Caprio|2005|l. 


4kpc, skewed toward larger values. Nonetheless, the 2kpc 
distance suggested if SGR 0418 indeed resides in the Perseus 
arm cannot firmly excluded. 

Finally, the best-fit is obtained for the following set of 
parameters: Nh ,22 ~ 0.24, tp = 96.8°, ^ = 7.9°, e = 0.08, 
n = 16961, Tpoie.oo = 0.38 keV and d = 2.27kpc, corre¬ 
sponding to a fit statistic xS/dof (prob.) = 58.5/38 (0.02). 
However, it is crucial to keep in mind that the full posterior 


distributions are the true representation of the acceptable 
parameter space given the model chosen and the data. 

Overall, this technique is superior to the iterative pro¬ 
cess of Section [4. l| since it includes all the effects of parame¬ 
ter variations on both the pulse profiles and spectra, which 
proved too difficult to perform iteratively. Nonetheless, the 
iterative analysis performed initially provided a validation 
of the simultaneous analysis via MCMC. 
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Figure 10. {left) Similar to Figure]^ but for a symmetry temperature distribution representing two antipodal spots. The solid line 
represents the probability p = 0.01. (right) Projected map of maximum p-values in a n—il> or n-J space. 


4.3 Can two antipodal spots be responsible for 
the observed pulse profile of SGR 0418? 

We present here an analysis similar to that of Section |4.1| 
but with a symmetric temperature distribution at the sur¬ 
face of the magnetar (Equation !^, i. e., with two antipodal 
spots. As can be seen in Figure ]l0] (left panel), the two- 
spots temperature distribution is able to accomodate the 
observed pulse profile for a narrow range of angles. One can 
note the additional symmetry (compared to the single spot 
case) since the two spots are identical. The two angles tp 
and ^ are constrained to be ^ 15° and ^ 165°, although 
the value of one angle strictly constrains the value of the 
other (and the two angles are symmetric with respect to ex¬ 
change). As with the single spot case, e 0.15, i.e. e is not 
constrained more than with the spectral analysis. 

However, we observed a strong correlation between the 
spots size and the angle geometry. Specifically, the larger 
the two spots, the more restricted the angles (see Figure [lOl 
right panel). This can easily be explained by the fact that 
two large spots (small n) can only reproduce the single- 
peaked pulse profile observed for angles tp and ^ with values 
~ 45° or ~ 135°. However, as the spots get smaller, they 
can accommodate the observed pulse fraction and profiles 
for wider ranges of angles. This results from the combined 
effect of beamed emission and light bending which allow the 
observer to see one spot appearing from behind the neutron 
star while the second spot is still visible. 

In this case of symmetric hot spots, the best fit to 
the pulsed profile is obtained for the parameters ('!/’,?) = 
(57°,30°), and e = 0.05 and n = 18000. As in Section [4.1[ 
we perform a consistency check using these parameters for a 
spectral fit. We obtain the best fit Tpoie.oo = 0.35±0.05 keV, 
Nh ,22 ~ 0.271q'o8i S' corresponding distance d = 1.51q 5 kpc, 
and a fit statistic x?/dof (prob.) = 1.26/34 (0.15). Figure]^ 
shows the best-fit pulse profile obtained with the symmetric 
temperature distribution at the surface of SGR 0418, along¬ 
side with the best-fit synthetic pulse profiles for the one-spot 
models (beaming and isotropic local emissions). 

Overall, we find that two identical small spots with 
a beamed local emission can produce high PFs and single 


peaked pulse profiles. Therefore, observed pulse profiles can 
be fitted equally well with a symmetric temperature distri¬ 
bution (two spots) as with asymmetric models (single spots), 
hence making the two situations equivalent. This result is 


similar to that derived for other sources (Shabaltas & Lai 
M^IStorch et al |2014| . 


5 INTERPRETATION OE OUR FINDINGS 
WITHIN THE CONTEXT OF HEATING 
MODELS 


The strongest constraint of our modelling is the presence of 
a high contrast in the temperature distribution on the sur¬ 
face of the star. The hottest point on the star, Tpoie.tx), is 
inferred to be at a value of about 0.35 keV. The tempera¬ 
ture then declines to an antipoidal minimum which has been 
constrained to be at least a factor of ~ 6 smaller than the 
maximum, but which can be much smaller (note that emis¬ 
sion from the coolest part of the star is not detectable in 
our observations). The precise rate of decline between the 
maximum and minimum values is not well determined by 
the current data. However, for the range of acceptable tem¬ 
perature profiles and distances, there is a region around the 
pole with an angular size of a few degrees for which the 
temperature hovers above ~ 0.3 keV. 

The inferred high temperature of the region dominating 
the emission is difficult to reconcile with the standard cool¬ 
ing model of NSs, even considering magneto-thermal evolu¬ 
tionary models with extra heating by magnetic field decay. 
The estimated age of the star, consistent with its luminosity 


and spin evolution, is about half a million years (Rea et al. 


20131. At that age, the temperatures of the hottest regions 


in a normal pulsar are expected to be below 0.1 keV, with 
some dependence on the equation of state, NS mass, su¬ 
perfluid gap, envelope model, etc. This result also holds for 
the available models with very s trong initial poloidal fields 
Br, ~ 10^® G (Vigano et al.|20]^l, or weak dipolar field with 


extremely strong initial toroidal fields, Bt ~ 10^® G (Gep- 
|pert fc Vigano|2014[ ), which is qualitatively more compatible 
with the timing properties and the outburst activity of SGR 
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0418. With some tuning of the microphysics parameters, the 
expected temperatnre of the hot spot at 0.5 Myr may be in¬ 
creased np to a 50% (say 0.15 keV). The mismatch between 
theory and observation conld be somewhat compensated by 
atmospheric effects, which may give a color correction by 
about a factor ~ 2 (e.g. Lloyd|[2C103 1, and hence could pos¬ 
sibly account for inferred temperatnres up to 0.3 keV. How¬ 
ever, the innermost region of the spot, of about a degree 
in size and temperatnre above 0.3 keV cannot be explained 
by residual cooling and internal heating in a 0.5 Myr old 
NS. We note that rather high temperatures accompanied by 
small inferred radii of the emitting regions are a common 
issue to most magnetars. In relatively young objects with 
strong toroidal fields and high A^h , magnetothermal simnla- 
tions show that km-sized hot spots can in fact be prodnced 
(Perna et al. 20131, bnt this is problematic in the case of 
SGR 0418 since it is much older than the rest of the sources, 
and the size of the hottest region is especially small. 

There are different possibilities to explain the anoma¬ 
lously high temperature. The first one is simply that the ob¬ 
ject has not reached its quiescence state yet; however, since 
the observations have showed a steady flux for more than 
a year, we believe that quiescence has indeed been reached. 
Therefore, unless something is missing from our theoretical 
understanding of magnetized NS cooling | |Pons et al.||200^ 
or from the physics of the NS crust and envelope, we con¬ 
clude that the hotspot must be maintained by an external 
heating source, attributed to energetic particles carried by 
magnetospheric currents and falling into the polar cap. The 
bombardment of these particles could account for the large, 
persistent temperature with the caveat of how to maintain 
stable current systems on a timescale of years. At present 
only analytical estimates of the bombardment energetics are 
available, without any detailed numerical simulation. Thus, 
deeper investigations are needed to conclude in favour of one 
or another option. 


6 CONCLUSIONS 

We have modelled the post-outburst surface emission of 
SGR 0418-1-5729. The low-luminosity emission, observed 
with XMM-Newton in 2012 and 2013 is consistent with being 
constant on this time-scale, hence indicating that the source 
has reached quiescence. Using a general-relativistic, phase- 
dependent thermal spectral model, we have fit the spectrum 
and pulse prohle of SGR 0418 to constrain the geometry of 
the system as well as the temperature distribution profile on 
the surface of the NS. 

This analysis was performed using two independent 
analyses: one approach which requires iteration between 
spectral analysis and pulse profile htting to constrain the 
parameters of the model; a second method using a Markov- 
Chain Monte Garlo approch that simultanesously htted the 
pulse prohle and spectra of SGR 0418. The two methods led 
to consistent results. 

We have found that SGR 0418 has a high temperature 
contrast on the surface, with differences between the maxi¬ 
mum and the minimum temperature of a factor of at least 
~ 6. Despite the single-peaked pulse prohle, the possibility 
of a symmetric temperature distribution at the surface of 
SGR 0418 (i.e., two antipodal spots) cannot however be ex¬ 


cluded. The small size of the spots, combined with radiation 
beaming and the absorption of soft X-rays, allow for a high 
PF single-peak pulse prohle observed for SGR 0418. Signif¬ 
icant constraints were also placed on the vie wing/emission 
geometry. While each angle can take any value between 0 
and 180°, we constrained ip and ^ in a correlated manner 
such that 65° < V' + ? < 125° or 235° < V' + C < 295° 
(with a mild dependence on the radiation beaming; for 
isotropic emission the requirement would be more constrain¬ 
ing: 104° < V’ + ? < 145° or 215° < V' + ^ < 256°). 

The inferred value of the NS surface temperature, ex¬ 
ceeding 0.3 keV in a region of about a few degrees in size, 
is very difhcult to explain in a source of about 0.5 Myr 
with standard cooling models, even accounting for possi¬ 
ble color corrections due to atmospheric processing of the 
emitted radiation, and varying the micro and macro physics 
parameters of the currently available cooling models. While 
a contribution from internal heating in the cooler region of 
the star cannot be excluded, we believe that the dominant, 
small hot spot must be maintained by the bombardament 
of energetic particles carried by magnetospheric currents. 
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